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ABSTRACT 


The  Thermodynamical  Ocean  Prediction  System  (TOPS)  Is  a 
general  and  flexible  software  framework  for  operational 
Implementation  of  upper  ocean  forecast  models  at  Fleet 
Numerical  Oceanography  Center  (FLENUMOCEANCEN),  Monterey, 
California.  It  was  developed  by  NORDA  Code  322  as  a  part  of 
the  Navy's  Automated  Environmental  Prediction  System  (AEPS). 
TOPS  Is  fully  Interfaced  with  the  FLENUMOCEANCEN  operational 
data  base  and  meets  all  FLENUMOCEANCEN  programming  standards 
for  operational  programs.  This  technical  note  discusses  the 
uses  for  TOPS,  provides  documentation  of  the  physics 
currently  represented  In  the  system,  Indicates  probable 
future  developments,  and  briefly  addresses  the  problem  of 
forecast  verification. 
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I.  INTRODUCTION 


Almost  universally,  the  upper  ocean  is  characterized  by  a  mixed-layer 
extending  from  the  surface  to  about  5-100  m  depth, in  which  temperature  and  salinity 
exhibit  only  small  changes  with  depth.  In  addition  to  having  important 
implications  for  long-range  weather  prediction,  climate  modeling,  fisheries 
operations,  ocean  thermal  energy  conversion  and  pollution  control,  the  depth  and 
stratification  of  this  layer  have  important  impacts  on  the  propagation  of 
underwater  sound.  Thus,  synoptic  knowledge  of  the  structure  of  the  oceanic  mixed 
layer  is  of  great  interest  to  the  Navy. 

The  mixed-layer  owes  its  high  degree  of  vertical  uniformity  to  mixing  caused 
by  turbulence  that  is  generated  by  shear  instabilities,  breaking  waves,  and  surface 
cooling.  A  dynamically  stable  water  mass  in  which  the  vertical  eddy  fluxes  are 
extremely  small  usually  exists  below  the  mixed-layer. 

During  periods  of  strong  wind  forcing  and/or  strong  surface  cooling,  the 
mixed-layer  tends  to  deepen  because  water  is  entrained  into  the  layer  from  below  by 
turbulent  mixing.  During  periods  of  relatively  weak  wind  forcing  and/or  strong 
surface  heating,  however,  the  source  of  turbulent  kinetic  energy  may  become  too 
weak  to  maintain  active  entrainment  at  the  base  of  the  mixed-layer,  causing  the 
layer  to  retreat  to  a  shallower  depth.  In  this  way,  the  thermal  structure  of  the 
mixed-layer  is  modified  substantially  on  time  scales  of  a  few  days  by  the  passage 
of  atmospheric  disturbances  (e.g.,  Elsberry  and  Camp,  1978;  Elsberry  and  Raney 
1978).  Hence,  the  temporal  variability  of  the  thermal  structure  of  the  upper  ocean 
is  much  larger  than  that  of  the  region  in  and  below  the  main  thermocline. 

Modeling  of  the  oceanic  mixed-layer  is  intrinsically  linked  to 
parameterization  of  turbulent  processes  and,  consequently,  is  a  difficult  problem. 
Over  the  past  decade,  however,  tremendous  interest  has  been  shown  in  upper  ocean 
modeling  and  much  progress  has  been  made.  For  example,  even  though  much  is  still 
unknown  about  the  fundamental  nature  of  turbulence,  field  data  have  shown  the  state 
of  the  mixed-layer  to  be  highly  predictable  with  a  variety  of  turbulence 
parameterization  models.  As  a  result,  it  is  now  reasonable  to  expand  the  realm  of 
mixed-layer  modeling  from  the  experimental  domain  to  the  operational  domain  by 
constructing  upper  ocean  forecast  systems  that  are  interfaced  with  operational  data 
bases.  The  Thermodynamical  Ocean  Prediction  System  (TOPS),  developed  for  Fleet 
Numerical  Oceanography  Center  (FLENUMOCEANCEN)  by  NORDA  Code  322,  is  the  first  such 
system.  TOPS  is  regarded  as  a  component  of  the  Navy's  steadily  evolving  Automated 
Environmental  Prediction  System  (AEPS). 

Several  important  uses  for  TOPS  are  anticipated.  For  example,  it  can  be  used 
to  improve  the  daily  Fleet  Numerical  Ocean  Thermal  Structure  Analysis  (which 
provides  initial  conditions  for  the  forecast  model)  by  producing  a  24-hour  forecast 
for  use  as  the  first-guess  field  (i.e.,  the  best  estimate  of  the  analyzed  field 
before  new  data  is  assimilated)  in  the  following  day's  analysis.  In  this  way,  the 
results  predicted  by  the  forecast  model  are  fed  back  into  the  analysis  on  a  daily 
basis.  This  should  Improve  the  analysis  everywhere  by  tending  to  make  it 
dynamically  consistent  with  the  atmospheric  forcing,  which  presumably,  is  known 
fairly  well.  The  improvement  should  be  especially  large  in  data-sparse  regions 
where,  in  the  absence  of  a  forecast  model ,  the  analyzed  thermal  field  is 
constrained  to  stay  very  near  climatology,  regardless  of  the  local  atmospheric 
forcing. 

In  addition  to  requiring  knowledge  of  the  Instantaneous  thermodynamical  state 
of  the  upper  ocean,  the  Navy  needs  to  have  the  capability  to  forecast  this  state 


over  periods  of  several  days,  as  this  ability  could  be  tactically  significant.  For 
example,  the  solar  heating  cycle  modifies  the  stratification  of  the  upper  ocean  in 
an  acoustically  important  way,  but  this  information  is  obviously  not  resolved  by 
the  daily  analysis.  Furthermore,  intense,  rapidly  moving  meteorological  distur¬ 
bances,  such  as  extratropical  cyclones  and  their  associated  warm  and  cold  fronts, 
can  substantially  modify  the  upper  ocean  over  vast  areas  during  a  period  of  72 
hours.  Given  proper  initial  conditions  and  a  reasonably  accurate  72-hour  weather 
prediction,  TOPS  can  provide  a  useful  72-hour  forecast  of  the  response  of  the  upper 
ocean  to  these  forcing  mechanisms. 

The  following  section  gives  a  detailed  description  of  the  physics  represented 
in  the  present  configuration  of  TOPS.  Section  III  summarizes  the  system  and  its 
uses,  gives  a  limited  view  of  probable  future  developments,  and  comments  briefly  on 
the  problem  of  forecast  verification. 

II.  PRESENT  CONFIGURATION  OF  TOPS 


Two  separate  forecast  models  are  available  in  the  present  configuration  of 
TOPS:  a  non-advective  or  "quasi-one-dimensional "  model  and  an  advective  or  "quasi 
three-dimensional"  model.  Both  models  use  the  same  parameterizations,  grid, 
initial  conditions,  and  upper  and  lower  boundary  conditions. 


A.  PROGNOSTIC  EQUATIONS  FOR  THE  NON-ADVECTIVE  MODEL  (HMLMS) 

In  the  non-advective  forecast  model  or  Hemispheric  Mixed  Layer  Model 
System  (HMLMS),  planetary  rotation  and  the  convergence  of  vertical  eddy  and 
radiative  fluxes  are  assumed  to  be  the  only  processes  controlling  the  dynamics  of 
the  upper  ocean.  Under  these  assumptions,  the  conservation  equations  for 
temperature,  salinity,  and  momentum  are 
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where  T  is  the  temperature,  S  the  salinity,  u  and  v  the  x-  and  y-components  of  the 
current  velocity  (x  and  y  relative  to  the  grid,  see  Section  II. E),  w  the  z- 
component  of  current  velocity,  F  the  downward  flux  of  solar  radiation,  D  a  damping 
coefficient,  v  a  diffusion  coefficient,  f  the  Coriolis  parameter,  t  the  time  and  z 
the  vertical  coordinate  (positive  upward).  Spatial  averages  over  a  region  defined 
by  a  cell  in  the  horizontal  mesh  are  denoted  by  (  )  and  primes  indicate  departure 

from  these  averages.  Thus,  for  example,  the  quantity  w'S'  represents  the  vertical 
eddy  (i.e.,  turbulent)  flux  of  salinity. 


The  terms  involving  the  damping  coefficient  D  in  (3)  and  (4)  represent 
the  drag  force  caused  by  the  stress  at  the  base  of  the  mixed-layer  associated  with 
the  propagation  of  internal  wave  energy  away  from  the  wind-forced  region  (e.g., 
Pollard  and  Millard,  1970).  As  discussed  by  Niiler  and  Kraus  (1977),  this  drag 
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force  can  contribute  to  the  relatively  fast  attenuation  of  inertial  oscillations 
observed  in  the  mixed-layer.  We  use  the  value  D  ■  0.1  day 1 ,  which  is  within  the 
range  of  estimates  for  this  quantity. 

The  terms  involving  v  in  (l)-(4)  represent  very  weak  "background"  eddy 
diffusion  that  exists  even  below  the  mixed-layer.  We  set  v  *  1  cm2  $-1  and  note 
that,  on  time  scales  less  than  one  month,  the  model  Is  not  sensitive  to  the  value 
of  ¥  ,  provided  that  it  is  not  zero. 

Equations  (l)-(4)  are  solved  on  a  three-dimensional  grid  In  order  to 
produce  a  three-dimensional  forecast  of  the  state  of  the  upper  ocean.  Since  there 
are  no  horizontal  processes  considered  In  this  formulation,  however,  it  is  best 
thought  of  as  a  "quasl-one-dimenslonal"  model. 

Over  most  of  the  world  ocean,  purely  one-dimensional  mixing  processes 
account  for  a  major  portion  of  the  upper  ocean's  response  to  strong  atmospheric 
forcing  (e.g..  Camp  and  Elsberry,  1978).  Hence,  even  though  this 
quasi-one-dimensional  formulation  Is  quite  simple,  it  should  be  capable  of 
representing  a  significant  part  of  the  upper  ocean's  short  time  scale  response, 
provided  that  suitable  parameterlzatlons  of  the  vertical  eddy  fluxes  are  used. 

B.  PROGNOSTIC  EQUATIONS  FOR  THE  ADVECTIVE  MODEL  (AHMLMS) 

In  the  advectlve  forecast  model  or  Advective  Hemispheric  Mixed  Layer 
Model  System  (AHMLMS),  horizontal  and  vertical  advection  and  horizontal  diffusion 
of  temperature  and  salinity  are  Included  in  addition  to  the  radiative  and  vertical 
mixing  processes  of  the  non-advective  model  (see  Section  II. A).  As  a  result,  the 
conservation  equations  for  temperature  and  salinity  become 
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where  ua,  va  and  wa  are  the  x-,  y-,  and  z-components  of  the  advection  current  and  A 
is  a  horizontal  eddy  diffusion  coefficient  (A  *  108  cm2  s"l).  Definitions  of  the 
remaining  symbols  can  be  found  either  in  Section  II. A  or  In  the  Appendix. 

Horizontal  pressure  gradients  and  horizontal  advection  and  diffusion  are 
neglected  In  the  momentum  equations  In  the  advective  model.  Therefore,  these 
equations  still  take  the  form  of  Equations  (3)  and  (4)  of  the  non-advective  model. 

Equations  (3)  and  (4)  can  be  used  to  determine  the  Ekman  component  of  the 
advection  current.  Since  there  are  no  pressure  gradient  terms  In  (3)  and  (4), 
however,  they  cannot  be  used  to  determine  the  geostrophlc  component.  Therefore, 
even  though  advection  and  diffusion  of  temperature  and  salinity  occur  in  three 
dimensions  In  this  formulation,  It  Is  not  a  true  three-dimensional  model.  It  Is, 
perhaps,  best  termed  a  "quasi -three-dimensional"  model. 


At  present,  the  geostrophlc  component  of  the  advectlon  current  is 
calculated  diagnostically  from  a  climatological  data  base  (see  Section  II. H. 2).  In 
future  applications,  the  geostrophlc  current  will  be  provided  by  a  hydrodynamical 
model  with  high  horizontal  resolution. 

As  mentioned  in  the  previous  section,  one-dimensional  mixing  and 
radiative  processes  dominate  the  dynamics  of  the  mixed-layer  on  time  scales  of  a 
few  days  over  most  of  the  world's  oceans.  On  longer  time  scales,  however, 
advective  processes  can  make  important  contributions  to  the  heat  and  salinity 
budgets  of  the  upper  ocean.  Therefore,  assuming  that  it  is  supplied  with  realistic 
advection  currents,  the  advective  model  should  perform  better  than  the 
non-advective  model  on  long  time  scales. 

The  advective  model  will  be  installed  on  the  CYBER  203  computer  at 
FLENUMOCEANCEN  in  1981  and  will  become  the  primary  operational  model  in  TOPS.  The 
non-advective  model  is  currently  operative  on  the  CYBER  175  at  Monterey.  Even 
after  the  advective  model  becomes  operational,  the  non-advective  version  will 
remain  active  by  serving  as  back-up  model  to  be  run  on  the  CYBER  175  in  the  event 
that  the  CYBER  203  is  down.  This  back-up  capability  will  help  guarantee  the 
operational  availability  of  the  product  fields  generated  by  TOPS. 

C.  PARAMETERIZATION  OF  TURBULENT  MIXING 

The  Level-2  turbulence  closure  theory  of  Mellor  and  Yamada  (1974)  is 
currently  used  to  parameterize  the  vertical  eddy  fluxes  of  temperature,  salinity 
and  momentum  in  both  the  advective  and  non-advective  models.  In  this  parameteri¬ 
zation,  the  fluxes  are  given  by 
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where  Kh  and  are  eddy  diffusion  coefficients,  t  is  the  turbulence  length  scale, 
q  is  the  square  root  of  twice  the  turbulent  kinetic  energy,  and  Sh  and  S^  are 
functions  of  the  gradient  Richardson  number  R1 ,  where 

Ri  =  pwJ»  .  (11) 

(ir  *  inn 

Here,  g  is  the  acceleration  of  gravity  and  p  is  the  mean-field  density  calculated 
from  T  and  S  according  to  the  equation  of  state  proposed  by  Freidrlch  and  Levltus 
(1972). 
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The  functional  forms  for  Sh  and  Sm  are  derived  by  Mel  lor  and  Yamada 
(1974)  and  require  specification  of  three  empirical  constants.  These  constants, 
however,  are  determined  once  and  for  all  from  neutrally  stratified  turbulent  flow 
data.  The  resulting  curves  for  Sh  and  are  shown  in  Figure  1.  Note  the  implied 
cutoff  of  turbulence  at  Ri  =  0.23. 

The  quantity  q  is  calculated  from  a  form  of  the  turbulent  kinetic  energy 
equation  that  expresses  a  local  balance  of  shear  production,  buoyancy  production, 
and  viscous  dissipation  of  turbulent  kinetic  energy 

£qSMl(ll)  +  (It)  1  *  A<,Sh(p^  IS")  ‘  15T=  0  (12 

where  ( 7 ) - ( 10 )  and  the  equation  of  state  have  been  used.  Hence,  the  basic 
assumption  of  this  turbulence  model  is  that  transports  of  turbulent  kinetic  energy 
can  be  neglected. 

Finally,  we  follow  Mel  lor  and  Durbin  (1975)  and  calculate  the  turbulence 
length  scale  from  the  ratio  of  the  first  to  the  zeroth  moment  of  the  turbulence 
field  /*o 


This  equation,  plus  (7)-(12),  closes  the  turbulence  parameterization.  This 
turbulence  model  has  been  applied  to  the  oceanic  mixed-layer  with  success  by  Mel  lor 
and  Durbin  (1975),  Martin  (1976),  Martin  and  Roberts  (1977),  Martin  and  Roberts 
(1978),  Clancy  (1979),  Martin  and  Thompson  (1980),  and  Warn-Varnas  et  al .  (1980). 

It  has  been  compared  to  higher-order  turbulence  closure  models  with  favorable 
results  by  Mel  lor  and  Yamada  (1974)  and  Warn-Varnas  and  Piacsek  (1979). 

D.  PARAMETERIZATION  OF  THE  EXTINCTION  OF  SOLAR  RADIATION 

Although  only  a  small  amount  of  the  incident  solar  radiation  penetrates 
below  the  upper  few  meters  of  the  sea  (about  70%  is  absorbed  in  the  upper  5  m, 
Jerlov,  1968),  the  penetration  of  solar  radiation  can  affect  mixed-layer  dynamics 
and  the  density  structure  of  the  upper  ocean  in  several  ways. 

1.  Penetration  of  solar  radiation  provides  a  means  of  warming  the  region 
below  the  mixed-layer  where  vertical  diffusion  of  heat  is  relatively  weak.  This 
process  is  especially  effective  when  the  mixed-layer  is  shallow. 

2.  Penetration  of  solar  radiation  allows  for  a  higher  level  of  turbulent 
kinetic  energy  within  the  mixed-layer.  Since  less  heat  has  to  be  mixed  downward, 
less  turbulent  kinetic  energy  is  absorbed  by  potential  energy  production.  This  can 
result  in  weaker  stability  at  the  base  of  the  mixed-layer  and  a  deeper  mixed-layer. 

3.  Penetration  of  solar  radiation  allows  convection  to  occur  between  the 
surface  and  the  compensation  depth  when  the  loss  of  heat  at  the  surface  is  less 
than  the  absorption  of  solar  radiation.  The  compensation  depth  is  defined  as  the 
depth  above  which  the  amount  of  solar  radiation  absorbed  by  the  sea  Is  equal  to  the 
loss  of  heat  from  the  surface.  Convection  can  occur  due  to  the  upward  flux  of  heat 
from  the  region  above  the  compensation  depth  to  the  surface. 

The  extinction  rate  of  solar  radiation  used  In  TOPS  is  from  Jerlov 
(1968).  Table  I  lists  the  percentage  of  solar  radiation  that  penetrates  to  certain 
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Figure  1.  The  stability  functions  S„  and  as  a  function  of  gradient 
Richardson  number  R^. 
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depths  for  several  types  of  water  in  the  open  ocean.  The  water  types  are 
classified  according  to  their  optical  clarity  from  Type  I  for  very  clear  ocean 
water  to  Type  III  for  fairly  turbid  ocean  water  such  as  right  be  found  in  a 
biologically  productive  region.  Figure  2  shows  the  worldwide  distribution  of  the 
various  water  types. 

The  divergence  of  the  solar  radiative  flux  in  TOPS  is  calculated  from 
il  =  p 

az  o  dz  (14) 

where  F0  is  the  total  solar  radiation  penetrating  the  surface  and  y  is  the  fraction 
of  F0  penetrating  to  a  depth  z.  The  quantity  v  is  interpolated  from  the  data  in 
Table  1. 

TOPS  currently  employs  an  extinction  profile  for  Type  IA  (ch-u'')  ocean 
water,  which  is  (according  to  Figure  2)  the  most  common  type  in  the  oper  ocean. 
However,  plans  are  to  amend  the  system  to  account  for  regional  and  perhaps  seasonal 
differences  in  water  type  according  to  Figure  2  and  more  recent  data. 

E.  GRID 

The  vertical  grid  used  in  both  the  advective  and  non-advecti ve  models 
consists  of  17  levels  between  the  surface  and  500  m  depth  and  is  shown  in  Figure  3. 
High  resolution  is  achieved  near  the  surface  in  order  to  properly  resolve  the 
mixed-layer,  and  every  fixed  level  in  the  FLENUMOCEANCEN  Expanded  Ocean  Thermal 
Structure  (EOTS)  analysis  is  represented  in  this  grid.  Note  that  the  vertical  eddy 
fluxes  and  wa  are  defined  at  depths  midway  between  those  for  which  temperature, 
salinity,  and  momentum  are  defined. 

For  the  horizontal  representation  in  both  models,  T,  S,  u  and  v  are 
defined  at  the  points  of  the  standard  FLENUMOCEANCEN  63  X  63  Northern  Hemisphere 
Polar  Stereographic  Grid  (see  Figure  4).  This  grid  is  true  at  60°N  where  the 
spacing  is  381  km. 

In  the  advective  model,  wa  is  also  defined  on  this  grid,  but  ua  and  va 
are  staggered  with  respect  to  these  points.  Figure  5  shows  the  basic  elements  of 
the  resulting  grid  system.  Note  that  in  both  model  formulations,  the  x-  and  y- 
directions  are  taken  relative  to  the  grid,  not  the  earth. 

F.  INITIAL  CONDITIONS 

The  initial  temperature  field  is  taken  from  the  daily  ocean  thermal 
structure  analysis  produced  by  the  FLENUMOCEANCEN  Operational  System.  TOPS  can 
accept  input  from  either  the  Expanded  Ocean  Thermal  Structure  (EOTS)  analysis 
system  or  the  older  Ocean  Thermal  Structure  ( OTS )  analysis  system. 

These  analysis  schemes  contain  no  explicitly  modeled  physics  and  are 
based  entirely  on  standard  Information  blending  concepts.  The  daily  Northern 
Hemisphere  data  set  input  to  these  systems  consists  of  about  200  XBT  observations, 
2000  bucket  sea-surface  temperature  observations,  and  20,000  satellite  sea-surface 
temperature  observations  from  TIROS-N.  Since  information  is  blended  vertically  as 
well  as  horizontally,  the  sea-surface  temperature  observations  contribute 
information  to  the  subsurface  thermal  analysis.  In  data-sparse  regions,  the 
analyzed  thermal  field  remains  very  near  a  state  determined  by  a  daily 
Interpolation  of  monthly  climatology.  This  monthly  climatology  is  based  on  an 


/ 


TABLE  1 

PERCENTAGE  OF  TOTAL  SURFACE 
IRRADIANCE  PENETRATING  TO  A  GIVEN 
DEPTH  FOR  A  SOLAR  ALTITUDE  OF  90° 
(from  Jerlov,  1968) 


Oceanic  Water 

Depth 

(m) 

I 

IA 

IB 

II 

III 

0 

100.0 

100.0 

100.0 

100.0 

100.0 

1 

44.5 

44.1 

42.9 

42.0 

39.4 

2 

38.5 

37.9 

36.0 

34.7 

30.3 

5 

30.2 

29.0 

25.8 

23.4 

16.8 

10 

22.2 

20.8 

16.9 

14.2 

7.6 

20 

25 

13.2 

11.1 

7.7 

4.2 

0.97 

50 

5.3 

3.3 

1.8 

0.70 

0.041 

75 

1.68 

0.95 

0.42 

0.124 

0.0018 

100 

0.53 

0,28 

0.10 

0.0228 

150 

0.056 

0.00080 

200 


0.0062 


Figure  4.  Standard  FLENUMOCEANCEN  63  x  63_Northern  Hemisphere  Polar 
Stereographic  Grid  on  which  T,  S,  u,  v,  and  wa  are  defined 


y 


Figure  5.  Staggered_gr1d  system  used  in  the  advective  model.  The  quantities 
T,  5,  u,  v,  and  wa  are  defined  at  the  points  denoted  •,  which  cor¬ 
respond  to  the  points  in  the  standard  FLENUMOCEANCEN  63  x  63  Northern 
Hemisphere  Polar  Stereographic  Grid.  The  quantity  ua  is  defined  at 
points  denoted  by  a.  and  va  is  defined  at  points  denoted  by  □.  The 
interpolated  climatological  density  field  and  the  geostrophic  stream 
function  are  defined  at  the  corners  of  the  dashed-line  boxes. 


12 


/ 


objective  analysis  of  approximately  600,000  bathythermograph  observations  made  in 
the  northern  hemisphere  prior  to  1974  (see  Weigle  and  Mendenhall,  1974). 

At  present,  the  first-guess  field  for  each  day's  analysis  is  generated 
from  the  previous  day's  analysis  by  a  forecast  of  persistence,  plus  adjustment 
toward  the  daily  interpolated  climatology.  As  mentioned  in  Section  I,  TOPS  will 
eventually  be  used  to  generate  the  first-guess  field,  and  thereby  bring  physics 
into  the  analysis  by  tending  to  make  it  dynamically  consistent  with  the  atmospheric 
forcing.  See  Holl  et  al .  (1979)  for  a  detailed  discussion  of  the  mathematical 
techniques  currently  used  in  the  analysis  schemes. 

The  initial  salinity  field  is  taken  from  a  daily  interpolation  of  monthly 
climatology.  Inclusion  of  salinity  in  the  forecast  models  is  necessary,  since  it 
makes  an  important  contribution  to  the  density  stratification  in  some  regions  and 
thereby  affects  the  vertical  turbulent  mixing.  Lack  of  a  synoptic  analysis  for 
salinity,  however,  can  present  a  problem  in  a  few  areas.  For  example,  anomalous 
(i.e.,  non-climatological )  salinity  distributions  can  sometimes  allow  strong 
temperature  inversions  to  exist  at  the  base  of  the  mixed-layer  in  high  latitudes. 

If  such  an  inversion  is  picked  up  by  the  daily  ocean  thermal  structure  analysis, 
then  an  unrealistic  convective  adjustment  will  occur  in  the  model  at  that  location. 
Use  of  temperature-salinity  (T-S)  relationships  is  not  an  adequate  solution  to  this 
problem,  since  these  relationships  cannot  be  applied  with  confidence  to  the  mixed- 
layer  in  these  anomalous  regions. 

To  overcome  the  problem,  an  option  exists  in  TOPS  which  allows  the 
salinity  to  be  adjusted  slightly  from  climatology  in  order  to  guarantee  that:  (1) 
the  initial  density  stratification  below  the  mixed-layer  does  not  fall  below  some 
user-specified  minimum  value,  and  (2)  the  initial  vertical  gradient  of  density  in 
the  mixed-layer  for  a  forecast  performed  on  day  N  is  given  by  the  24  hour  forecast 
vertical  density  gradient  in  the  mixed-layer  calculated  on  day  N-l.  To  begi n  a 
sequence  of  daily  forecast  runs,  the  salinity  Is  adjusted  to  make  the  initial 
density  stratification  neutral  in  the  mixed-layer. 

The  initial  momentum  field  is  also  provided  by  carrying  information 
forward  in  time  from  day  to  day.  To  begin  a  sequence  of  daily  forecast  runs,  the 
initial  momentum  profile  is  set  to  zero  below  the  base  of  the  mixed-layer  and  taken 
to  vary  linearly  in  the  mixed-layer  such  that  the  mixed-layer-averaged  flow  is 
equal  to  the  vertically  averaged  Ekman  drift.  In  successive  forecast  runs, 
however,  the  initial  momentum  field  is  given  by  the  24  hour  forecast  momentum  field 
calculated  by  the  previous  day's  run. 

The  turbulence  length  scale  4  (see  Section  II.C)  is  initialized  in  a 
similar  manner.  To  begin  a  sequence  of  forecast  runs,  4  is  set  to  2.5  m.  In 
subsequent  forecast  runs,  4  is  provided  by  the  value  produced  by  the  24  hour 
forecast  of  the  previous  day's  run. 

Note  that  with  the  initial  conditions  supplied  in  the  manner  described 
above,  the  turbulent  kinetic  energy,  the  vertical  eddy  diffusion  coefficients,  the 
vertical  eddy  flux  of  momentum,  and  the  vertical  eddy  flux  of  density  (to  the 
approximation  that  a  linear  equation  of  state  holds)  are  conserved  in  the  process 
of  updating  the  model -predicted  temperature  field  with  the  dally  ocean  thermal 
structure  analysis  (see  Section  II.C).  This  treatment  tends  to  make  the  Initial 
state  consistent  with  both  the  atmospheric  forcing  and  the  dynamics  of  the 
turbulence  parameterization  model,  which  lessens  the  "initialization  shock" 
associated  with  bringing  the  temperature  observations  into  the  model  each  day  via 
the  daily  analysis,  and  restarting  the  forecast. 


G.  BOUNDARY  CONDITIONS 

The  upper  boundary  conditions  for  the  temperature,  salinity,  and  momentum 
conservation  equations  are  provided  by  surface  fluxes  that  are  forecast,  out  to  a 
period  of  72  hours  on  the  same  grid  used  by  TOPS,  by  operational  FLENUMOCEANCEN 
atmospheric  models.  Thus, 
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where  B0  is  the  surface  infrared  radiative  flux,  H0  is  the  surface  sensible  heat 
flux,  LQ0  is  the  surface  latent  heat  flux,  P0  is  the  surface  precipitation  rate,  S0 
is  the  surface  salinity,  andrx  andrY  are  the  components  of  the  surface  wind  stress. 
Definitions  for  the  remaining  symbols  can  be  found  in  either  earlier  sections  or  in 
the  Appendix. 

The  components  of  the  surface  wind  stress  are  computed  at  6  hour 
intervals  from 

/  =  »aCD  U{U2  +  V2)*  ,  (19, 
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where  U  and  V  are  the  x-  and  y-components  of  the  6-hourly  available  wind  velocity 
vector  at  a  reference  level  above  the  sea  surface,  pa  is  a  reference  density  for 
air,  and  Cq  is  a  constant  drag  coefficient.  The  resulting  components  of  the  wind 
stress  and  the  quantities  (B0  +  H0  +  IQ«)  and  Q0,  which  are  also  derived  from  the 
FLENUMOCEANCEN  fields  at  6  hour  intervals,  are  then  Interpolated  to  each  time  step 
(0.5  hour)  of  the  forecast  models.  At  present,  only  a  simple  linear  Interpolation 
scheme  is  used.  However,  a  higher  order  Interpolation  scheme  will  be  adopted  in 
the  future,  if  it  proves  warranted. 

The  surface  precipitation  rate  P0  is  derived  from  a  predicted  field  of  12 
hour  accumulated  precipitation.  It  is  simply  assumed  constant  during  each  12  hour 
period  of  the  forecast. 

The  surface  solar  radiation  flux  F0  (which  provides  the  upper  boundary 
condition  for  the  radiational  heating  calculation)  is  available  from  the 
atmospheric  model  only  at  6  hour  Intervals.  It  is  Interpolated  to  each  time  step 
of  the  oceanic  models  according  to 

F0  =  I  cosa  (21) 
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where  a  is  the  time-varying  local  zenith  angle  of  the  sun  and  I  is  a  linear 
function  of  time  chosen  such  that  the  resulting  curve  for  F0  passes  through  all  of 
the  6-hourly  available  values.  This  treatment  is  necessary  to  allow  adequate 
representation  of  the  solar  flux. 

The  studies  of  Johnson  (1977),  Paul  us  (1978)  and  Elsberry  et  al .  (1979) 
have  indicated  that  the  wind  stress  and  heat  flux  fields  derived  from  the 
FLENUMOCEANCEN  operational  system  do  provide  fairly  realistic  synoptic  time  scale 
forcing  parameters  for  ocean  prediction.  However,  problem  areas,  where  the 
integrated  surface  heat  flux  is  not  consistent  with  the  change  in  the  heat  content 
of  the  upper  ocean  derived  from  observations,  have  been  identified  (see  Elsberry  et 
al.,  1979).  This  difficulty  seems  to  be  most  severe  south  of  latitude  30°  and  is 
probably  related  to  the  fact  that  the  present  atmospheric  model  is  only  hemispheric 
in  coverage.  Implementation  of  the  Naval  Operational  Global  Atmospheric  Prediction 
System  (NOGAPS)  will  probably  improve  the  surface  flux  fields  in  this  region 
considerably. 

The  lower  boundary  conditions  for  the  conservation  equations  are  provided 
by  holding  the  temperature,  salinity,  and  momentum  constant  at  the  lower  boundary 
of  the  model  during  each  forecast  run.  The  temperature  and  salinity  are  defined 
there  by  a  linear  extrapolation  of  the  initial  conditions  downward  below  400  m. 

The  momentum  field  at  the  bottom  boundary  is  always  set  to  zero. 

Because  there  are  no  horizontal  exchanges  of  heat  and  salinity  in  the 
non-advective  model,  no  lateral  boundary  conditions  are  required  in  that 
formulation.  In  the  advectivejnodel _,  the  normal  component  of  the  advection  current 
and  the  normal  derivatives  of  T  and  S  are  taken  to  be  zero  at  land-sea  boundaries. 
Thus,  no  advection  or  diffusion  of  heat  or  salinity  is  allowed  across  these 
boundaries.  In  addition,  the  normal  derivatives  of  T  and  S  on  the  outer  boundary 
surrounding  the  forecast  domain  (i.e.,  one-half  grid  space  outside  of  the  63  X  63 
grid)  are  assumed  to  be  zero,  which  implies  no  diffusion  of  heat  and  salinity  into 
or  out  of  the  domain.  Finally,  the  corner  ua  values  (i.e.,  ua  (1,1),  ua  (1,63),  ua 
(64,1),  ua  (64,63))  are  set  equal  to  the  nearest  interior  ua  values  and  the 
vertical  component  of  the  advection  current  is  assumed  to  be  zero  along  the  open 
boundary  of  the  63  X  63  grid.  This  allows  the  horizontal  components  of  the 
advection  current  to  be  calculated  from  continuity  along  the  outer  boundary. 

H.  CALCULATION  OF  THE  ADVECTION  CURRENT 

The  ocean  current  used  to  advect  the  temperature  and  salinity  fields  in 


the  forecast  model 

is  given  by 

ua  -  ue  +  u§  , 

and 

va  *  ve  +  vg  » 

wa  =  we  , 

(22) 

where  ue,  ve,  and  we  are  the  x-,  y-,  and  z-components  of  the  wind-driven  Ekman 
circulation  and  uA  and  vA  are  the  x-  and  y-components  of  a  nondlvergent  geostrophlc 
velocity  field  determined  from  a  climatological  data  base. 
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1 .  Ekman  Component  of  the  Advection  Current 

Only  the  steady  component  of  the  wind-driven  Ekman  current  which  is 
in  balance  with  the  wind  stress  and  the  mixed-layer  depth  is  used  for  advection. 

The  time-dependent  Ekman  current  calculated  as  part  of  the  mixed-layer  model 
(Equations  3  and  4)  is  not  used  in  order  to  filter  out  inertial  oscillations. 
Because  of  their  periodic  nature,  inertial  oscillations  are  not  very  effective  in 
advecting  the  density  field  over  distances  on  the  scale  of  the  model  grid  which  has 
a  resolution  of  200  to  400  km.  For  example,  the  diameter  of  a  20  cm  s_1  inertial 
current  circle  at  30°N  is  only  17  km.  Filtering  the  inertial  fluctuations  from  the 
Ekman  current  allows  the  use  of  a  large  time  step  for  advection  and  a  consequent 
saving  of  computer  time.  The  advection  terms  are  updated  every  6  hours  during  the 
forecast,  which  is  often  enough  to  resolve  changes  in  the  wind  field,  but  not  often 
enough  to  resolve  inertial  oscillations  without  introducing  noise  into  the 
horizontal  and  vertical  advection  fields. 


The  wind-driven  Ekman  current  is  calculated  using  the  equations 
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where  the  symbols  have  the  same  meanings  defined  previously. 


The  boundary  conditions  for  these  equations  are  that  the  surface 
stress  is  equal  to  the  wind  stress  and  that  the  stress  at  the  base  of  the  mixed 
layer  (z  *  -h)  is  zero: 


(25) 

(26) 

(27) 

(28) 


The  zero  stress  condition  at  the  base  of  the  mixed-layer  is  a  good  approximation, 
since  turbulent  mixing  in  this  region  is  almost  completely  suppressed  by  the 
stratification  and  the  stress  is  quite  small  there  relative  to  the  mean  stress 
within  the  mixed-layer. 


The  surface  wind  stress  and  the  mixed-layer  depth  are  required  to 
define  the  boundary  conditions.  Forecast  values  of  the  wind  stress  are  provided  by 
FLENUMOCEANCEN  atmospheric  models  (as  discussed  in  Section  II. G)  and  the 
mixed-layer  depth  is  obtained  from  TOPS  itself. 

As  in  the  mixed-layer  equations,  the  Mellor-Yamada  Level- 2  turbulence 
closure  scheme,  described  in  Section  II. C,  is  used  to  parameterize  the 
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turbulent  momentum  fluxes  w'u'  and  w V .  However,  for  the  calculation  of  the 
steady  Ekman  advection  current,  the  mixed-layer  is  taken  to  be  completely  mixed 
(i.e.,  unstratified  with  a  Richardson  number  of  zero).  Hence,  the  momentum  fluxes 
are  a  function  only  of  the  vertical  shear  of  the  horizontal  velocities  ue  and  ve. 

Equations  (23)  and  (24)  are  solved  successively  for  each  column  of 
the  grid  at  the  T-S  points  (see  Fig.  5).  The  Ekman  velocities  ue  and  ve  are  then 
horizontally  interpolated  to  the  points  where  the  advection  velocities  ua  and  va 
are  defined.  The  vertical  grid  used  is  as  described  in  Section  II. E.  As  an 
initial  guess  to  the  Ekman  velocity  profiles,  (23)  and  (24)  are  solved  directly  by 
setting  the  vertical  eddy  diffusion  coefficient  for  momentum  equal  to  50  cm2  s-l. 
The  time  dependent  forms  of  (23)  and  (24)  (see  Equations  (3)  and  (4))  are  then 
integrated  with  a  one-hour  time  step  to  convergence,  with  the  eddy  coefficients 
being  updated  each  time  step  using  the  Mellor-Yamada  turbulence  parameterization. 


The  vertical  motion  that  results  from  the  horizontal  divergence  of 
the  Ekman  current  field  is  calculated  by  integrating  the  continuity  equation 


9X 


3y 


5 

3Z 


=  0 


from  the  surface  to  a  depth  z  giving 
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(30) 


The  vertical  velocity  at  the  surface  has  been  taken  to  be  zero.  In  actual 
calculations,  we  is  defined  at  the  same  points  as  wa  (see  Section  II. E). 

The  vertical  motion  just  below  the  mixed-layer  due  to  the  divergence 
of  the  surface  Ekman  current  is  sometimes  referred  to  as  the  Ekman  pumping  or  Ekman 
suction  velocity.  By  substituting  the  steady  Ekman  equations  into  the  continuity 
equation  (29)  and  ignoring  the  drag  terms  (which  are  small)  we  can  obtain 
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In  Equation  (31),  the  Ekman  pumping  velocity  we  (-h)  depends  only  on  the  wind 
stress  curl  and  to  a  lesser  extent  on  the  latitudinal  variation  of  the  Coriolis 
parameter  f.  It  is  notable  that  the  rate  of  Ekman  pumping  does  not  depend  upon  the 
particular  parameterization  used  for  turbulent  mixing  or  the  mixed-layer  depth,  and 
is  therefore  independent  of  the  mixed-layer  model  itself. 

2.  Geostrophlc  Component  of  the  Advection  Current 

The  geostrophic  currents  used  for  advection  are  calculated  from  the 
FLENUMOCEANCEN  climatological  temperature  and  salinity  fields.  These  fields 
are  monthly  averages  and  are  available  on  the  FLENUMOCEANCEN  63  X  63  Northern 
Hemisphere  Polar  Stereographic  Grid  used  by  TOPS  at  standardized  depths  from  the 
surface  to  5000  m. 

Consideration  was  given  to  using  the  daily  ocean  thermal  structure 
analysis  to  calculate  the  geostrophlc  advection  currents.  However,  a  trial 
calculation  In  which  the  climatological  temperature  field  In  the  upper  400  m  was 
replaced  by  the  current  OTS  analysis  resulted  in  a  significant  weakening  of  certain 
currents  compared  to  the  climatological  calculation.  Although  the  mid-latitude 


currents  such  as  the  Gulf  Stream  and  Kuroshlo  looked  similar,  the  equatorial 
current  system  In  the  Pacific  was  severely  diminished.  The  reason  for  this 
difference  is  unclear,  but  it  suggests  that  the  sparse  data  input  to  the  analysis 
in  this  region  is  insufficient  to  resolve  the  horizontal  density  gradients. 
Because  of  this  problem,  the  climatological  temperature  field  is  presently  being 
used  to  calculate  the  geostrophic  component  of  the  advection  current. 


The  procedure  for  calculating  the  geostrophic  component  of  the 
advection  current  is:  (1)  the  density  field  is  calculated  from  the  appropriate 
monthly  temperature  and  salinity  fields,  (2)  geostrophic  currents  are  calculated 
from  the  density  field  using  the  thermal  wind  relations  and  a  reference  level  of  no 
motion,  and  (3)  a  stream  function-vorticity  equation  is  solved  to  eliminate  the 
horizontal  divergence  of  the  geostrophic  current  and  to  satisfy  the  lateral 
boundary  condition  of  no  flow  across  land-sea  boundaries. 


The  density  field  is  calculated  from  the  monthly  climatological 
temperature  and  salinity  fields  using  a  polynomial  formulation  of  the  equation  of 
state  for  seawater  developed  by  Friedrich  and  Levitus  (1972).  The  densities  are 
calculated  at  the  T-S  points  of  the  grid  where  the  climatological  temperature  and 
salinity  fields  are  defined,  and  then  horizontally  interpolated  to  the  corner 
points  of  the  dashed-line  boxes  of  Figure  5.  This  allows  the  geostrophic 
velocities  to  be  directly  calculated  at  the  points  where  the  advection  velocities 
ua  and  va  are  located. 


The  thermal  wind  equations  are  used  to  calculate  the  geostrophic 
current  from  the  density  field.  These  relations  (given  by  Pond  and  Pickard,  1978) 
are 


(32) 


and 
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The  use  of  the  above  equations  allows  only  the  calculation  of  relative  currents. 

In  order  to  calculate  absolute  currents,  the  absolute  velocity  at  some  depth  must 
be  specified.  An  assumption  frequently  used  is  that,  at  a  certain  depth,  the 
currents  become  small  and  can  be  taken  to  be  zero.  Such  an  assumption,  when  used 
to  estimate  current  transports  over  a  deep  water  column,  can  lead  to  significant 
errors,  since  deep  currents  can  be  large  enough  to  yield  appreciable  transports 
when  integrated  over  a  large  depth.  Since  deep  currents  tend  to  be  much  weaker  than 
surface  currents,  however,  the  use  of  a  depth  of  no  motion  to  estinate  upper-ocean 
currents  generally  yields  reasonably  small  errors. 

The  reference  level  or  level  of  no  motion  that  was  selected  for 
calculating  geostrophic  currents  from  the  climatological  density  fields  is  a 
function  of  latitude  and  is  based  on  a  number  of  trial  calculations  using  reference 
levels  between  400  and  2500  m.  Below  the  equator  (l.e.,  in  the  "corners"  of  the  63 
X  63  northern  hemisphere  grid)  the  deep  density  climatology  Is  extremely  poor  and  a 
reference  level  of  500  m  Is  used.  The  use  of  a  deeper  reference  level  In  this 
region  generates  spurious  currents.  Above  30°N  a  reference  level  of  1250  m  Is 
used.  The  strongest  currents  In  this  region,  the  Gulf  Stream  and  the  Kuroshlo, 
have  depth  scales  of  about  1000  m.  These  currents  are  not  much  enhanced  by  using  a 
deeper  reference  level,  but  are  significantly  reduced  when  the  reference  depth  is 
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decreased.  Between  0°N  and  30°N  the  reference  depth  was  arbitrarily  taken  to  be  a 
linear  function  of  latitude  to  provide  a  smooth  transition  of  the  level  of  no 
motion  between  these  two  latitudes.  It  was  found  that  the  equatorial  currents 
calculated  from  the  climatological  density  fields  between  0°N  and  30°N  are 
relatively  shallow  and  are  not  sensitive  to  the  depth  of  the  reference  level  as 
long  as  it  is  greater  than  about  400  m.  The  use  of  a  reference  level  for  the 
calculation  of  geostrophic  currents  that  increases  with  latitude  is  consistent  with 
the  findings  of  previous  investigators.  For  example,  Defant  (1941)  found  the 
optimum  reference  level  to  be  about  600  m  near  the  equator,  increasing  to  1250-2000 
m  at  40°N  in  the  Atlantic. 

A  problem  occurs  when  applying  the  thermal  wind  equations  near  the 
equator  where  f  approaches  zero.  Although  the  geostrophic  approximation  is 
considered  to  be  fairly  good  to  within  one  or  two  degrees  of  the  equator,  small 
errors  in  the  density  field  near  the  equator  can  generate  large  spurious  currents. 
For  this  reason,  the  value  of  f  used  for  calculating  the  geostrophic  current  in  the 
region  between  5°N  and  5°S  is  taken  to  be 


f  ■  f  (at  5°N)  for  0°  <  latitude  <  5°N 


and 


(34) 


f  =  f  (at  5°S)  for  5°S  <  latitude  <  0° 

This  treatment  yields  a  realistic  looking  equatorial  current  system. 

It  is  desirable  to  eliminate  the  horizontal  divergence  of  the 
geostrophic  component  of  the  advection  current.  Since  calculating  divergence 
involves  taking  a  small  difference  between  large  terms,  small  errors  in  the 
horizontal  current  (due  to  noise  in  the  density  field)  tend  to  cause  large  errors 
in  the  vertical  motion  field,  which  is  calculated  from  the  divergence  of  the 
horizontal  geostrophic  current.  Since  the  geostrophic  current  is  calculated  from 
monthly  climatology,  erroneous  features  in  the  vertical  motion  field  will  persist 
for  at  least  a  month  until  this  field  is  updated.  It  TOPS  is  used  for  long 
forecasts,  or  if  the  model  forecast  is  fed  back  into  the  daily  ocean  thermal 
structure  analysis,  these  persistent  errors  in  the  vertical  advection  field  could 
have  cumulative  effects  that  would  generate  significant  errors  in  the  predicted  and 
analyzed  density  fields. 


In  order  to  eliminate  the  horizontal  divergence  of  the  geostrophic 
velocity  field,  a  stream  function-vorticity  equation  is  solved  at  each  level.  This 
equation  is  given  by 


lit  +  lit  =  ?!a .  !!!a  , 

9x2  ay2  3y 


(35) 


where  the  right-hand-side  of  the  above  equation  is  the  vertical  component  of 
vorticity  of  the  divergent  geostrophic  velocity  field  and  t  is  the  stream  function. 
A  horizontally  non-divergent  velocity  field  is  obtained  from  the  stream  function 
using  the  standard  definitions 


(36) 


and 


(37) 


Here,  the  velocity  field  defined  from  the  stream  function  is  designated  by  an 
asterisk  to  differentiate  it  from  the  divergent  geostrophic  field  that  appears  on 
the  left-hand-side  of  (32)  and  (33). 

The  stream  function  is  defined  at  the  corners  of  the  dashed-line 
boxes  shown  in  Figure  5.  On  the  land-sea  boundaries,  the  stream  function  is  set  to 
zero  so  that  there  will  be  no  flow  across  these  boundaries.  The  single  exception 
is  that  the  stream  function  on  the  Cuba-Haiti  island  group  is  defined  to  be  (1  - 
depth/1000  m)  x  5X10^  m2  s_l  to  provide  a  transport  of  25  Sverdrups  (25  X  10°  m3 
s‘l)  through  the  Yucatan  and  Florida  .Straits.  Since  the  Florida  Straits  span  only 
one  grid  interval,  this  is  the  simplest  means  of  getting  a  reasonable  flow  througn 
this  region. 


Boundary  conditions  along  the  open  boundaries  of  the  grid  which  span 
the  South  Atlantic,  the  South  Pacific,  and  the  Indian  Ocean  (see  Figure  4)  are 
defined  by  normalizing  the  geostrophic  transport  across  each  of  these  boundaries  to 
zero  by  means  of  an  additive  constant.  The  normalization  of  the  transport  across 
the  open  boundaries  has  the  effect  of  distorting  the  flow  near  some  parts  of  the 
boundary.  This  distortion  is  most  severe  along  the  South  Atlantic  boundary  because 
the  near-surface  geostrophic  flow  across  this  boundary  is  almost  entirely 
northward.  However,  the  distortion  diminishes  significantly  within  10  to  20°  of 
the  boundary. 


Equation  (35)  is  solved  for  the  stream  function  at  each  level  using 
successive-over-relaxation.  Because  the  model  grid  is  relatively  coarse,  the 
solution  converges  sufficiently  in  just  a  few  iterations. 

III.  SUMMARY  AND  OUTLOOK 


The  Thermodynamical  Ocean  Prediction  System  (TOPS)  is  a  general  and  flexible 
software  framework  for  operational  implementation  of  upper  ocean  forecast  models  at 
Fleet  Numerical  Oceanography  Center  (FLENUMOCEANCEN) .  It  was  developed  by  NORDA 
Code  322  as  a  part  of  the  Navy's  Automated  Environmental  Prediction  System  (AEPS). 
TOPS  is  fully  interfaced  with  the  FLENUMOCEANCEN  operational  data  base  and  meets 
all  FLENUMOCEANCEN  programming  standards  for  operational  programs. 

The  horizontal  grid  used  by  TOPS  is  the  standard  FLENUMOCEANCEN  63  X  63 
Northern  Hemisphere  Polar  Stereographic  Grid.  The  vertical  grid  consists  of  17 
levels  between  the  surface  and  500  m  depth  with  stretching  employed  to  retain  high 
resolution  in  the  upper  100  m  for  proper  treatment  of  the  mixed-layer. 

TOPS  is  initialized  by  the  daily  FLENUMOCEANCEN  analysis  of  ocean  thermal 
structure  and  a  daily  interpolation  of  monthly  climatological  salinity  fields.  The 
ocean  predictions  are  driven  out  to  a  forecast  time  of  72  hours  by  fluxes  of  heat, 
moisture  and  momentum  at  the  sea  surface  supplied  by  operational  FLENUMOCEANCEN 
atmospheric  models. 

Two  forecast  models  are  available  in  the  system.  In  the  non-advective  model, 
the  time  rate  of  change  of  temperature  and  salinity  is  due  only  to  the  convergence 
of  vertical  eddy  and  radiative  fluxes.  Thus,  even  though  the  model  produces  a 
three-dimensional  forecast  of  the  upper  ocean,  it  contains  only  one-dimensional 
physical  processes  and  is,  therefore,  termed  a  "quasi-one-dimensional"  model.  In 
the  advective  model ,  horizontal  and  vertical  advection  and  horizontal  diffusion  of 
temperature  and  salinity  are  included  in  addition  to  the  radiative  and  vertical 
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mixing  processes  of  the  non-advective  model.  The  horizontal  pressure  gradient  term 
is  still  not  included  in  the  equations  of  motion;  therefore,  the  model  cannot 
be  used  to  determine  the  geostrophic  component  of  the  flow.  Thus,  even  though 
advection  and  diffusion  of  temperature  and  salinity  occur  in  three  dimensions  in 
this  formulation,  it  is  not  a  true  three-dimensional  model  and  is  therefore  termed 
a  "quasi-three-dimensional"  model. 

In  the  present  configuration  of  the  advective  model,  the  geostrophic 
component  of  the  advection  current  is  calculated  diagnostically  from  a 
climatological  data  base.  The  Ekman  component  of  the  advection  current  is 
determined  from  the  predicted  local  values  of  mixed  layer  depth  and  surface  wind 
stress. 


In  both  models,  the  Level -2  turbulence  closure  theory  of  Mel  lor  and  Yamada 
(1974)  is  used  to  parameterize  the  vertical  eddy  fluxes  of  heat,  salinity,  and 
momentum.  The  basic  assumption  of  this  parameterization  is  that  the  turbulent 
kinetic  energy  budget  at  each  level  in  the  mixed-layer  consists  of  a  balance 
between  shear  production,  buoyancy  production  and  viscous  dissipation. 

The  advective  model  will  be  installed  on  the  CYBER  203  computer  at 
FLENUMOCEANCEN  in  1981  and  will  become  the  primary  operational  model  in  TOPS.  The 
non-advective  model  is  currently  operative  on  the  CYBER  175  at  Monterey,  and  even 
after  the  advective  model  becomes  operational,  will  remain  active  by  serving  as  a 
back-up  model  to  be  run  on  the  CYBER  175  in  the  event  that  the  CYBER  203  is  down. 
This  back-up  capability  is  important  in  an  operational  context  and  will  help 
guarantee  the  daily  availability  of  the  product  fields  generated  by  TOPS. 

Since  the  thermodynamical  structure  of  the  upper  ocean  has  important 
implications  for  many  activities,  there  are  multiple  uses  for  TOPS.  From  the 
Navy's  viewpoint,  however,  the  uses  for  TOPS  are  mainly  twofold:  (1)  improve  the 
daily  analysis  of  upper  ocean  thermal  structure  by  generating  24-hour  forecasts  for 
use  as  first-guess  fields  in  the  analysis,  and  (2)  provide  real-time  72-hour 
forecasts  of  tactically  significant  changes  in  the  thermal  structure  of  the 
mixed- layer. 

From  a  software  standpoint,  TOPS  has  been  designed  in  a  highly  modular 
fashion.  This  will  allow  the  system  to  undergo  continual  evolution  with  relative 
ease.  In  future  applications,  for  example,  TOPS  will  be  coupled  with  a 
hydrodynamical  model  (presently  under  development  by  NORDA  Code  322)  that  will 
supply  geostrophic  advection  currents  that  are  more  dynamically  consistent  than 
those  now  used.  Since  proper  treatment  of  hydrodynamical  processes  requires 
adequate  representation  of  mesoscale  eddies  and,  hence,  extremely  high  horizontal 
resolution,  the  hydrodynamical  model  will  necessarily  have  coarse  vertical 
resolution.  Thus,  although  the  hydrodynamical  component  of  this  coupled 
hybrid-system  will  have  coarse  vertical  resolution,  the  thermodynamical  component 
will  retain  fine  vertical  resolution  in  order  to  allow  proper  representation  of 
thermodynamical  processes  and  provide  meaningful  input  to  acoustic  models. 
Furthermore,  zoomed  (i.e.,  limited  area,  fine  horizontal  resolution)  versions  of 
the  system  may  be  produced  to  complement  the  zoomed  versions  of  the  daily  ocean 
thermal  structure  analysis  that  are  available.  In  addition,  more  sophisticated 
turbulence  parameterization  models  may  eventually  be  utilized. 

At  all  stages  in  the  development  of  operational  forecast  systems,  forecast 
verification  must  play  a  significant  role.  Although  the  turbulence 
parameterization  scheme  currently  used  in  the  system  (which  forms  the  heart  of  the 
forecast  models)  has  been  tested  successfully  in  a  number  of  one-dimensional 


studies  by  several  authors  (see  Section  II. C),  it  is  still  quite  necessary  to  test 
the  overall  forecast  system.  This  is  true  mainly  because  (1)  TOPS  depends  on  the 
FLENUMOCEANCEN  operational  data  base  for  initial  conditions  and  surface  forcing; 
consequently,  its  value  is  closely  linked  to  the  quality  of  this  data,  (2) 
performing  in  operational  mode,  the  turbulence  model  will  be  subjected  to  a  range 
of  conditions  wider  than  what  was  considered  in  the  one-dimensional  studies,  and 
(3)  in  the  advective  model,  three-dimensional  processes  are  also  included. 

During  the  developmental  stage,  informal  and  limited  test  and  evaluation 
studies  were  conducted  on  TOPS  by  NORDA  personnel  with  very  encouraging  results. 
Formal  and  extensive  test  and  evaluation  programs  on  the  present  and  future 
versions  of  TOPS  will  be  carried  out  primarily  by  FLENUMOCEANCEN  personnel  and 
represent  important  steps  before  the  various  models  can  achieve  full  operational 
status. 

Quantitative  verification  of  a  large-scale  ocean  forecast  poses  special 
problems  because  of  the  nonuniform  and  constantly  changing  spatial  distribution  of 
the  limited  number  of  daily-available  XBT  and  sea-surface  temperature  observations. 
Meaningful  forecast  verification  can  be  accomplished  only  in  data-rich  regions, 
since  the  model  forecast  will  probably  be  more  realistic  than  the  subsequent 
analysis  in  data-sparse  regions.  This  is  the  case  because  the  forecast  model  will 
respond  to  the  atmospheric  forcing  in  a  dynamically  consistent  way  everywhere, 
while  the  analysis  will  simply  remain  near  the  climatological  state  in  data-sparse 
regions,  regardless  of  the  local  atmospheric  forcing.  Thus,  when  these  test  and 
evaluation  programs  are  performed,  the  distribution  of  the  various  observations 
must  be  monitored  on  a  day-by-day  basis.  This  will  allow  subregions  of  the  grid 
that  have  adequate  data  coverage  at  both  the  beginning  and  end  of  a  forecast  period 
to  be  identified.  Model  forecasts  in  these  subregions  can  then  be  quantitatively 
compared  to  forecasts  of  persistence  and  climatology,  with  the  daily  analysis  at 
the  end  of  the  forecast  period  providing  the  verification  data  (the  model  forecast 
would  not  be  used  to  generate  the  first-guess  field  for  the  analysis  in  these 
comparisons).  Finally,  when  evaluating  an  ocean  prediction  system  in  this  way,  the 
relative  contribution  of  the  three  primary  sources  of  apparent  forecast  error,  (1) 
improper  initial  conditions,  (2)  inaccurate  atmospheric  forcing,  (3)  imperfect 
representation  of  oceanic  dynamics  and  (4)  inexact  verification  data,  must  be 
sorted  out. 


IV.  REFERENCES 


Camp,  N.T.  and  R.L.  Elsberry  (1978).  Oceanic  Thermal  Response  to  Strong  Atmospher¬ 
ic  Forcing  II.  The  Role  of  One-Oimensional  Processes.  J.  Phys.  Oceanog.,  8, 
215-224. 

Clancy,  R.M.  (1979).  A  Model  of  Diurnal  Variability  of  the  Ocean-Atmosphere  System 
in  the  Undisturbed  Trade  Wind  Regime.  Tech.  Report  SAI-79-807-WA,  Science 
Applications,  Inc.,  8400  Westpark  Drive,  McLean,  VA,  July,  114  p. 

Defant,  A.,  (1941).  Die  Absolute  Topographic  des  Physikalishen  Meersniveaus  and 
der  DruckflHchen,  Sowie  die  Wasserbewegungen  im  Atlantischen  Ozean.  Meteor- 
Werk,  6,  No.  2,  Liefg.  5,  191-260. 

Elsberry,  R.L.  and  N.T.  Camp  (1978).  Oceanic  Thermal  Response  to  Strong  Atmospher¬ 
ic  Forcing  I.  Characteristics  of  Forcing  Events.  J.  Phys.  Oceanog.,  8, 
206-214. 

Elsberry,  R.L.  and  S.D.  Raney  (1978).  Sea  Surface  Temperature  Response  to  Varia¬ 
tions  in  Atmospheric  Wind  Forcing.  J.  Phys.  Oceanog.,  8,  881-887. 

Elsberry,  R.L.,  P.C.  Gallacher,  and  R.W.  Garwood  (1979).  One-Dimensional  Model 
Predictions  of  Ocean  Temperature  Anomalies  During  Fall  1976.  Tech  Report  NPS 
63-79-003,  Naval  Postgraduate  School,  Monterey,  CA,  August,  30  p. 

Freidrich,  H.,  and  S.  Levitus  (1972).  An  Approximation  to  the  Equation  of  State 
for  Sea  Water,  Suitable  for  Numerical  Ocean  Models.  J.  Phys.  Oceanog.,  2, 
514-517. 


Holi ,  M.M.,  M.J.  Cuming,  and  B.R.  Mendenhall  (1979).  The  Expanded  Ocean  Thermal- 
Structure  Analysis  System:  A  Development  Based  on  the  Fields  by  Information 
Blending  Methodology.  Tech.  Report  M-241,  Meteorology  International 
Incorporated,  205  Montecito  Avenue,  Monterey,  CA,  July,  216  p. 

Jerlov,  N.G.  (1968).  Optical  Oceanography.  New  York,  Elsevier,  352  p. 

Johnson,  W.F.  (1977).  Upper  Ocean  Thermal  Structure  Forecast  Evaluation  of  a 

Model  Using  Synoptic  Data.  M.S.  Thesis,  Naval  Post  Graduate  School,  Monterey, 
CA,  47  p. 

Martin,  P.J.  (1976).  A  Comparison  of  Three  Diffusion  Models  of  the  Upper  Mixed 
Layer  of  the  Ocean.  NRL  Memorandum  Report  3399,  Naval  Research  Laboratory, 
Washington,  DC,  62  p. 

Martin,  P.J.  and  G.O.  Roberts  (1977).  An  Estimate  of  the  Impact  of  OTEC  Operation 
on  the  Vertical  Distribution  of  Heat  in  the  Gulf  of  Mexico.  Proc.  of  4th  OTEC 
Symposium,  P.  IV-26,  G.  loup.  (ed  1). 

Martin,  P.J.  and  G.O.  Roberts  (1978).  Peak  Current  Profile  Estimates  at  Selected 
Sites  for  Ocean  Thermal  Energy  Conversion.  SAI  Tech.  Report,  Science 
Applications,  Inc.,  8400  Westpark  Drive,  McLean,  VA. 

Martin,  P.J.  and  J.D.  Thompson  (1980).  Formulation  and  Testing  of  a  Layer- 
Compatible  Upper  Ocean  Mixed-Layer  Model.  Paper  submitted  to  Journal  of 
Physical  Oceanography. 


w 


Mel  lor,  G.L.  and  T,  Yamada  (1974).  A  Hierarchy  of  Turbulence  Closure  Models  for 
Planetary  Boundary  Layers.  J.  Atmos.  Sci.,  31,  1791-1806. 

Mel  lor,  G.L.  and  P.A.  Durbin  (1975).  The  Structure  and  Dynamics  of  the  Ocean  Sur¬ 
face  Mixed  Layer.  J.  Phys.  Oceanog.,  5,  718-725. 

Niiler,  P.P.  and  E.B.  Kraus  (1977).  One-Dimensional  Models  of  the  Upper  Ocean. 
Chapt.  10,  Modelling  and  Prediction  of  the  Upper  Layers  of  the  Ocean.  New 
York,  Pergamon  Press,  325  p.  (Pergamon  Marine  Series,  v.  1). 

Paulus,  R.A.  (1978).  Salinity  Effects  in  an  Oceanic  Mixed-Layer  Model.  M.S. 
Thesis,  Naval  Post  Graduate  School,  Monterey,  CA,  /7  p. 

Pollard,  R.T.  and  R.C.  Millard  (1970).  Comparison  Between  Observed  and  Simulated 
Wind-Generated  Inertial  Oscillations.  Deep  Sea  Res.,  17,  813-821. 

Pond,  S.  and  G.L.  Pickard  (1978).  Introductory  Dynamic  Oceanography.  New  York, 
Pergamon  Press,  241  p. 

Warn-Varnas,  A.C.  and  S.A.  Piacsek  (1979).  An  Investigation  of  the  Importance  of 
Third-Order  Correlations  and  Choice  of  Length  Scale  in  Mixed  Layer  Modelling. 
Geophys.  Astrophys.  Fluid  Dyn.,  13,  225-243. 

Warn-Varnas,  A.C.,  G.M.  Dawson,  P.J.  Martin,  and  M.  Miyake  (1980).  Forecasts  and 
Studies  of  the  Oceanic  Mixed  Layer  During  the  MILE  Experiment.  Paper  Submitted 
to  Journal  of  Geophysical  and  Astrophysical  Fluid  Dynamics. 

Weigle,  W.  F.  and  B.  R.  Mendenhall  (1974).  Climatology  of  the  Upper  Thermal- 
Structure  of  the  Seas.  Tech.  Report  M-196,  Meteorology  International, 
Incorporated,  205  Montecito  Avenue,  Monterey,  CA,  79  p. 


24 


APPENDIX 


SYMBOL 

A 

Bo 

c 

Cd 

D  * 
F 

Fo 

f 

9 

h 

Ho 

I 

Kh 

km 

l 

L 

P0 

q 

Qo 

Ri 

S 

So 

SH 

Sm 


LIST  OF  SYMBOLS 

DEFINITION 

Horizontal  eddy  diffusion  coefficient. 

Upward  infrared  radiation  flux  at  sea  surface. 

Specific  heat  of  seawater. 

Drag  coefficient  for  surface  wind  stress  calculation. 

Damping  coefficient  for  inertial  oscillations. 

Downward  flux  of  solar  radiation. 

Downward  flux  of  solar  radiation  at  sea  surface. 

Coriolis  parameter. 

Acceleration  of  gravity. 

Thickness  of  mixed  layer. 

Upward  sensible  heat  flux  at  sea  surface. 

Downward  flux  of  solar  radiation  at  sea  surface  scaled  by  the 
cosine  of  the  local  zenith  angle. 

Vertical  eddy  diffusion  coefficient  for  heat  and  salinity. 
Vertical  eddy  diffusion  coefficient  for  momentum. 

Turbulence  length  scale. 

Latent  heat  of  evaporation  for  water. 

Surface  precipitation  rate. 

The  square  root  of  twice  the  turbulent  kinetic  energy. 

Upward  evaporative  flux  of  moisture  at  sea  surface. 

Gradient  Richardson  number. 

Salinity. 

Salinity  at  sea  surface. 

Stability  function  for  vertical  eddy  fluxes  of  heat  and  salinity. 
Stability  function  for  vertical  eddy  flux  of  momentum. 


t 


Time. 


T 

u 

ue 

u9 

ug 

ua 

U 

v 

ve 

v9 

vg 

va 

V 

w 

we 

wa 

x 

y 

z 

r ) 
O 

a 

y 


Temperature. 

x-component  of  current  velocity. 

x-component  of  steady  Ekman  part  of  advection  current. 

x-component  of  geostrophic  current. 

x-component  of  divergence-free  geostrophic  part  of  advection 
current. 

x-component  of  advection  current. 

x-component  of  wind  velocity  at  reference  height  above  sea 
surface. 

y-component  of  current  velocity. 

y-component  of  steady  Ekman  part  of  advection  current, 
y-component  of  geostrophic  current. 

y-component  of  divergence- free  geostrophic  part  of  advection 
current. 

y-component  of  advection  current. 

y-component  of  wind  velocity  at  reference  height  above  sea 
surface. 

z-component  of  current  velocity. 

Ekman-induced  z-component  of  current  velocity. 

z-component  of  advection  current. 

Grid-referenced  horizontal  coordinate  (see  Fig.  4). 

Grid- referenced  horizontal  coordinate  (see  Fig.  4). 

Vertical  coordinate,  positive  upward  from  sea  surface. 

Spatial  average  at  constant  depth  taken  over  region  defined  by 
grid  cell  of  horizontal  grid. 

Departure  from  above-defined  spatial  average. 

Local  zenith  angle  of  the  sun. 

Fraction  of  surface  flux  of  solar  radiation  pentrating  to  a  depth 
z. 

Background  vertical  eddy  diffusion  coefficient. 


Density  of  seawater. 

Reference  density  for  air. 

Reference  density  for  water. 

Stream  function. 

x-component  of  surface  wind  stress  vector 
y-component  of  surface  wind  stress  vector 
Surface  wind  stress  vector. 
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